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1 Introduction 

The purpose of this paper. This paper provides a method for taking advantage of continuous 
symmetries in solving Lie group invariant optimization problems for cost functions that are defined 
on k^^-oideT tangent spaces of Lie groups. The type of application we have in mind is, for example, 
the interpolation and comparison of a series of images in longitudinal studies in a biomedical setting. 
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Previous work on the geometric theory of Lagrangian reduction by symmetry on first-order 
tangent spaces of Lie groups provides a convenient departure point that is generahzed here to allow 
for invariant variational problems formulated on higher-order tangent spaces of Lie groups. It turns 
out that this generalization may be accomplished as a series of adaptations of previous advances 
in Euler-Poincare theory, placed into the context of higher-order tangent spaces. Extension of 
the basic theory presented here to allow for actions of Lie groups on Riemannian manifolds should 
have several interesting applications, particularly in image registration, but perhaps elsewhere, too. 
Actions of Lie groups on Riemannian manifolds will be investigated in a subsequent treatment. 
Two important references for the present work are [HMR98] for the basic Euler-Poincare theory 
and [CMROl] for the bundle setting of geometric mechanics. 

1.1 Previous work on geometric splines for trajectory planning and in- 
terpolation 

The topics treated here fit into a class of problems in control theory called trajectory planning 
and interpolation by variational curves. These problems arise in numerous applications in which 
velocities, accelerations, and sometimes higher-order derivatives of the interpolation path need to 
be optimized simultaneously. Trajectory planning using variational curves in Lie groups acting on 
Riemannian manifolds has been discussed extensively in the literature. For example, trajectory 
planning for rigid body motion involves interpolation on either the orthogonal group 5*0(3) of 
rotations in M^, or the semidirect-product group SE{3) ~ S'0(3)(S)]R^ of three-dimensional rota- 
tions and translations in Euclidean space. Trajectory planning problems have historically found 
great utility with applications, for example, in aeronautics, robotics, biomechanics, and air traffic 
control. 

Investigations of the trajectory planning problem motivated the introduction in [GK85] and 
[NHP89] of a class of variational curves called Riemannian cubics. Riemannian cubics and their 
recent higher order generalizations are reviewed in [PopOT] and [MSKIO], to which we refer for 
extensive references and historical discussions. The latter work addresses the interpolation by 
variational curves that generalizes the classical least squares problem to Riemannian manifolds. 
This generalization is also based on the formulation of higher-order variational problems, whose 
solutions are smooth curves minimizing the L^-norm of the covariant derivative of order A; > 1, 
that fit a given data set of points at given times. These solutions are called fc^'^-order geometric 
splines, or geometric fc-splines. This approach was initiated in [NHP89] for the construction of 
smoothing splines with k = 2 for the Lie group SO (3) and then generalized to higher order in 
[CSC95]. The following result, noted in the first of these papers and then discussed more generally 
in the second one, was another source of motivation for the present work. 

Proposition 1.1 ([NHP89]). 

The equation for a 2^"^ -order geometric spline for a bi-invariant metric on S0{3) may be written 
as a dynamical equation for a time- dependent vector f2(t) G M"^ using the vector cross product 



for all t in a certain interval [0,T]. 

Solutions of the more general version of equation (1.1) expressed in [CS95] for 2"'^-order geo- 
metric splines on Lie groups in terms of the Lie algebra commutator are called 'Lie quadratics' in 
[Noa03, Noa04, Noa06]. 
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As we said, understanding the intriguing result in Proposition 1.1 from the viewpoint of Lie 
group-invariant higher-order variational principles was one of the motivations for the present work. 
Its general version is proved again below as equation (3.21) in Section 3 by using the Euler-Poincare 
methods of [HMROn] for higher-order variational principles that are invariant under the action of a 
Lie group. The directness and simplicity of the present proof of the general version of Proposition 

1.1 compared with other proofs available in the literature encouraged us to continue investigating 
the application of Lie group-invariant k^^-oideT variational principles for geometric fc-splines. It 
turns out that higher-order Euler-Poincare theory is the perfect tool for studying geometric k- 
splines. 

The Euler-Poincare theory for first-order invariant variational principles focuses on the study of 
geodesies on Lie groups, which turns out to be the fundamental basis for both ideal fluid dynamics 
and modern large-deformation image registration. For reviews and references to earlier work on 
first-order invariant variational principles, see [HMR98] for ideal fluids and [YoulO] for large- 
deformation image registration. The present paper begins by extending these earlier results for 
geodesies governed by first-order variational principles that are invariant under a Lie group, so as to 
include dependence on higher-order tangent spaces of the group (i.e., higher-order time derivatives 
of curves on the group). This extension is precisely what is needed in designing geometric /c-splines 
for trajectory planning problems on Lie groups. The essential strategy in making this extension is 
the application of reduction by symmetry to the Lagrangian before taking variations, as introduced 
in [HMR98] for continuum dynamics. The equivalence of the result of Lagrangian reduction by 
symmetry with the results in the literature for Riemannian cubics and k^'^-ordei geometric splines 
is shown in Section 3, Proposition 3.3. 

This previous work has created the potential for many possible applications. In this paper, 
we shall concentrate on the application of these ideas in template matching for Computational 
Anatomy {CA). Although we do not perform explicit image matching here, we demonstrate the 
higher-order approach to template matching in the finite dimensional case by interpolating a 
sequence of points on the sphere S"^, using S0{3) as the Lie group of transformations. 

1.2 Main content of the paper 

The main content of the paper is outlined as follows: 

Section 2 discusses the geometric setting for the present investigation of extensions of group- 
invariant variational principles to higher order. In particular. Section 2 summarizes the 
definition of higher order tangent bundles and connection-like structures defined on them, 
mainly by adapting the treatment in [CMROl] for the geometric formulation of Lagrangian 
reduction. 

Section 3 explains the quotient map for higher-order Lagrangian reduction by symmetry and uses 
it to derive the basic A;*^-order Euler-Poincare equations. This extends to higher-order the 
Euler-Poincare equations derived in [HMR98]. The k'^^-oidei Euler-Poincare equations are 
then applied to derive the equations for geometric /c-splines on a Lie group. After these 
preliminary developments, there follows a sequence of adaptations of previous advances in 
Euler-Poincare theory to higher-order tangent spaces. 

Section 4 extends the Clebsch-Pontryagin approach of [GBR 10] to develop the k^^-oidei Euler- 
Poincare equations for potential applications in optimal control. This extension highlights 
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the role of coadjoint motion for cotangent-lift momentum maps. 

Section 5 addresses theoretical and numerical results for our main motivation, longitudinal data 
interpolation. That is, interpolation through a sequence of data points. After a brief account 
of the previous work done in Computational Anatomy (CA), we derive the equations that 
generalize the equations for geodesic template matching [BGBHRIO] to the case of higher- 
order cost functionals and sequences of several data points. We recover in particular the 
higher-order Euler-Poincare equations. For a particular choice of cost functionals one can 
therefore think of the higher-order template matching approach as template matching by ge- 
ometric k-splines. We discuss the gain in smoothness afforded by the higher-order approach, 
then we provide a qualitative discussion of two Lagrangians that are of interest for appli- 
cations in CA. Finally, we close the section by demonstrating the higher-order approach to 
template matching in the finite dimensional case by interpolating a sequence of points on 
the sphere S"^, using 5*0(3) as the Lie group of transformations. This yields the template- 
matching analog of the NHP equation of [NHP89] in (1.1). The results are shown as curves 
on the sphere in Figures 5.2. A sample figure is shown below to explain the type of results 
we obtain. 




Fig. 1.1: First order vs. second order template matching results interpolating a sequence of evenly time-separated 
points on the sphere, using a bi-invariant metric on the rotation group SO{3). The colors show the local speed 
along the curves on the spheres (white smaller, red larger). The motion slows as the curve tightens. 

Section 6 extends to fe'^^-order tangents the metamorphosis approach of [HTY09] for image regis- 
tration and the optimization dynamics introduced in [GBHRIO]. 

Section 7 addresses Hamiltonian and Hamilton-Ostrogradsky formulations of the higher-order 
Euler-Poincare theory. The Hamilton-Ostrogradsky formulation results in a compound Pois- 
son bracket comprising a sum of canonical and Lie-Poisson brackets. 

Section 8 discusses the outlook for future research and other potential applications of the present 
approach. These include the formulation of higher-order Lie group invariant variational 
principles that include both curves on Lie groups and the actions of Lie groups on smooth 
manifolds, and the formulation of a /cth-order brachistochrone problem. 
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This paper represents only the beginning of our work in this direction. The extensions to higher 
order discussed here demonstrate the unity and versatihty of the geometric approach. We hope 
these methods will be a source of inspiration for future analysis and applications of Lie group 
reduction of higher-order invariant variational problems. 



2 Geometric setting 

We shall begin by reviewing the definition of higher order tangent bundles and the connection- 
like structures defined on them. For more details and explanations of the geometric setting for 
higher-order variational principles see [CMROl]. 

2.1 k^'^-order tangent bundles 

The fc*''-order tangent bundle Tq'' : T^^^Q Q is defined as the set of equivalence classes of 
curves in Q under the equivalence relation that identifies two given curves qi{t),i = 1,2, if 
?i(0) = ^2(0) = go and in any local chart we have qf\o) = q2\o), for / = 1, 2, . . . , /c, where g*^'^ 
denotes the derivative of order /. The equivalence class of the curve q{t) at go ^ Q is denoted [gjgo''- 
The projection 

t^q'^-.TC'^Q^Q is given by r^^^Mj?) = go- 

It is clear that T(°)Q = Q, T^^^Q = TQ, and that, for < / < A;, there is a well defined fiber 
bundle structure 

^a,.) ^ ^ given by ^ ([g]£)) = [g]£. 

Apart from the cases where A; = and k = 1, the bundles T^^^Q are not vector bundles. The 
bundle T^'^^Q is often denoted Q, and is called the second order bundle. 

Remark 2.1. We note that T^'^^Q = Jq(R,Q) consists of fc-jets of curves from M to Q based at 
G M, as defined, for example, in [Bou71, §12.1.2]. 

A smooth map f : M ^ N induces a map 

TWf , T^^)M ^ T^AT given by T^/ ([g]^) := [/ o g]^';)^). (2.1) 

In particular, a group action $ : G x Q — ?■ Q naturally lifts to a group action 

$(^) : G X T(^)Q ^ T('=)Q given by $f ([g]i:)) := T^*, ([gU^'^) = [-^-p ° 'Zll'U • (2.2) 

This action endows T^'^^Q with a principal G-bundle structure. The quotient {T^^^Q) /G is a fiber 
bundle over the base Q/G. The class of the element [gjgo'' in the quotient (T'^'^^Q) /G is denoted 



(fc) 

90 



J G 



The case of a Lie group. The fc^'^-order tangent bundle T^^^G of a Lie group G carries a natural 
Lie group structure: if [g]^go ■, and [^j^^^ are classes of curves g and h in G, define [5']go''[/i]|^Q'' 



[S'^lgoL' '^^^ algebra TeT^^^G of T^^^G can be naturally identified, as a vector space, with 
{k + 1)0 (that is, the direct sum of + 1 copies of q) which, therefore, carries a unique Lie algebra 
structure such that this identification becomes a Lie algebra isomorphism. 
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2.2 k^'^-order Euler-Lagrange equations 

Consider a Lagrangian L : T^'^^Q — > M, L = L {q,q,q, Then a curve q : [to,ti] —^-Qisa. 

critical curve of the action 

J[q]= r L{q{t),q{t),-;<l'^'\t))dt (2.3) 

J to 

among all curves q{t) E Q whose first (k — l) derivatives are fixed at the endpoints: q^^\ti), i = 0,1, 
j = 0, k — 1, if and only if q{t) is a solution of the k^^-oidei Euler-Lagrange equations 

i:(-l)^'^#^=0- (2-4) 

^ dp dq(3) ^ ' 

The corresponding variational principle is Hamilton's principle, 

5 r L{q{t),q{t),-;q^'\t))dt = 0. 

Jto 

In the 5-notation, an infinitesimal variation of the curve q{t) is denoted by 6q{t) and defined by 
the variational derivative, 

d 



q(t,e), (2.5) 



e=0 



where q{t,0) = q{t) for all t for which the curve is defined and ^{ti,e) = q^^\ti), for all e, 
j = 0, 1, . . . , A; - 1, i = 0, 1. Thus 5q^^\to) = = 5q^^\tx) for j = 0, k - 1. Note that the local 

notation L (g, q, g'^'^^) used above can be intrinsically written as L f [q]^''^ 



Examples: Riemannian cubic polynomials and generalizations. As originally introduced 
in [NHP89], Riemannian cubic polynomials generalize Euclidean splines to Riemannian manifolds. 
Let {Q, 7) be a Riemannian manifold and be the covariant derivative along curves associated 
with the Levi-Civita connection V for the metric 7. The Riemannian cubic polynomials are defined 
as minimizers of the functional in (2.3) for the Lagrangian L : T(2)g ^ M defined by 

H<l,<i,Q) ■=l%(^J,^j] ■ (2.6) 



2 Dt ^ 

This Lagrangian is well-defined on the second-order tangent bundle since, in coordinates 

D 



Dt 



q=q' + ^{qW, (2.7) 



where 0^ij{q))i,j,k are the Christoffel symbols at point q of the metric 7 in the given basis. These 
Riemannian cubic polynomials have been generalized to the so-called elastic splines through the 
following class of Lagrangians 

1 D \ , , , 

Lr{q, g, q) ■■= ^Iq ( j^q, j^q j + y7<?(^, q), (2.8) 
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L, (g, q, gW) := -7, ^^g, y-^g , (2.9) 



where r is a real constant. Another extension are the higher-order Riemannian sphnes, or geometric 
A;-sphnes, where 

for k > 2. As for the Riemannian cubic sphnes, Lk is weh-defined on T^'^^Q. Denoting by R 
the curvature tensor defined as R{X,Y)Z = V^VyX — VyVx^ — ^[x,y]Z, the Euler-Lagrange 
equation for elastic splines (fc = 2) reads 

+ R f ^g(t), m] m = r^^m, (2.10) 



Dt^^'^ ' -^Dt^' ''^'^J^' ' Dt 

as proven in [NHP89]. For the higher-order Lagrangians L^, the Euler-Lagrange equations read 
[CSC95] 



i=2 



These various Lagrangians can be used to interpolate between given configurations on T^^^Q. 
The choice of Lagrangian will depend on the application one has in mind. For instance, the 
following interpolation problem was addressed in [HB04a] and was motivated by applications in 
space-based interferometric imaging. 



Interpolation problem. Given N+1 points qi G Q, i = 0, and tangent vectors Vj G Tg.Q, 
j = 0,N, minimize 

among curves t q(t) G Q that are on [to,t]^], smooth on [tj,tj_|.i], to < ti < . . . < t^, and 
subject to the interpolation constraints 

q{U) = Qh for all i = 1, . . . ,N - 1 

and the boundary conditions 

Q{to) = go, qito) = Vo, and q^t^) = q^, q{tN) = Vn- 

In the context of a group action and invariant Lagrangians, we refer the reader to Section 5 for an 
example of higher-order interpolation particularly relevant for Computational Anatomy. 



2.3 Quotient space and reduced Lagrangian 

When one deals with a Lagrangian L : T^^^Q — t- M that is invariant with respect to the lift 
^^'^^ : G X T^'^^Q — i- T^'^^Q of a group action $ : G x Q — t- then the invariance can be exploited 
to define a new function called the reduced Lagrangian on the quotient space (T^'^^Q^ /G. We 
review this procedure here. Since this paper mainly deals with the case where Q = G, we begin 
by describing this special case. 
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Let G be a Lie group and h E G. The right-, respectively left-actions by h on G, 

Rh : G ^ G, g ^ gh, and G ^ G, g ^ kg, 

can be naturally lifted to actions on the A;*^-order tangent bundle T^'^^G (see (2.2)). We will denote 
these lifted actions by concatenation, as in 



: T('=)G ^ T('^)G, [gW ^ Ri'H[9& =■■ [9]l% and 



Consider a Lagrangian L : T^^^G R that is right-, or left-invariant, i.e., invariant with respect 
to the lifted right-, or left-actions of G on itself. For any [g]l'l^ G T^'^'^G we then get 



(2.13) 



respectively. The restriction LI (fc) of the Lagrangian to the fc*''-order tangent space at the identity 
e therefore fully specifies the Lagrangian L. Moreover, there are natural identifications ak '■ 
T^^^G kQ given by 



or 



ock{[gt^) ■.= [m 



d 

dt 
di 



■k-l 



t=0 



g{t)-'g{t),---, 



df"-^ 



t=0 



t=0 



dt^-^ 



t=0 



(2.14) 



(2.15) 



respectively, where t (?(t) is an arbitrary representative of [g]e'^- 
The reduced Lagrangian £ : A;g — )■ M is then defined as 

i := L|^(fe)^ o a^^, 



(2.16) 



where one uses the choice for that is appropriate, namely (2.14) for a right-invariant Lagrangian 
L and (2.15) for a left-invariant Lagrangian L. Let 1 1— )■ g{t) G G be a curve on the Lie group. For 
every t this curve defines an element in T^^-^G, namely 



[g]^g^t) '■~ [^laW' '^here h is the curve r 1—^ /i(r) := g{t + r). 



(2.17) 



Note that for the case = 1 we write, as usual, g{t) := [g]^g(ty The following lemma is a direct 
consequence of the definitions: 



Lemma 2.2. Let t 1— )■ g(t) be a curve in G and L : T'^^^G 
grangian. Then the following equation holds for any time to, 



a right-, or left-invariant La- 

(2.18) 



L{[g]lt))=^{ato),iito),...,^^'-'\to)), 

where ^ := gg~^ , or ^ := g~^g respectively. 

This last equation will play a key role in the higher-order Euler-Poincare reduction discussed 
in the next section. 
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3 Higher-order Euler-Poincare reduction 

In this section we derive the basic A;*^-order Euler-Poincare equations by reducing the variational 
principle associated to the Euler-Lagrange equations on T^'^^Q. The equations adopt a factorized 
form, in which the Euler-Poincare operator at /c = 1 is applied to the Euler-Lagrange operation 
acting on the reduced Lagrangian ■ ■ ■ ,^^''~^^) : fcg — ?■ M at the given order, k. We then 

apply the fc*''-order Euler-Poincare equations to derive the equations for geometric fc-splines. 

3.1 Quotient map, variations and k^'^-order Euler-Poincare equations 

Let L : T^'^^G ^ M be a right-, or left-invariant Lagrangian. Recall from §2.2 that the Euler- 
Lagrange equations are equivalent to the following variational problem: 

For given hi ^ G and [/i] G Tji'^'^G, i = 1,2, find a critical curve of the functional 

J{9] = L {[g]%) dt 

among all curves g :t & [^1,^2] ^ g{t) € G satisfying the endpoint condition 

[9]%)^ = [h]t'\ ^ = 1,2. (3.1) 

The time derivatives of up to order k — 1 are therefore fixed at the endpoints, i.e., = 

j = 0, . . . , k — 1, are automatically verified. Let g : t g{t) G G be a curve and (e, t) 1— j- geif) G G 

a variation of g respecting (3.1). We recall from Lemma 2.2 that, for any e and any tg, 

^([^4!(to))=^fc(^o),...,er^Hio)), (3.2) 

where := geQj^ ■, or := g^^ge respectively for the right-, or left-invariant Lagrangian L. The 
variation 5^ induced by the variation 5g is given by 

K = VT[i.r,], (3.3) 

where rj := {Sg)g~^, or rj := g~^{Sg), respectively. It follows from the endpoint conditions (3.1) 
that 7]{ti) = ffiti) = ... = 7]^''-^\ti) = and therefore 6^{ti) = ... = dt^sliti) = 0, for i = 1,2. 
We are now ready to compute the variation of J': 
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were we used the vanishing endpoint conditions 6^{ti) = . . . = d^~^5^{ti) = and ri{ti) = 0, for 
i = 1,2, when integrating by parts. Therefore, the stationarity condition 6 J' = imphes the 
-order Euler-Poincare equation. 



{dt±^di)J2{-iydi— = o. (3.4) 

j=0 



Formula (3.4) takes the following forms for various choices of A; = 1, 2, 3: 
If A; = 1: 

{d, ± ad*) - = 0, 

lfk = 2: 



If A; = 3: 



(9,±ad*)(--9,^)=0, (3.5) 



The first of these is the usual Euler-Poincare equation. The others adopt a factorized form in 
which the Euler-Poincare operator [dt ± ad^ is applied to the Euler-Lagrange operation on the 
reduced Lagrangian ■■■) at the given order. 

The results obtained above are summarized in the following theorem. 

Theorem 3.1 (/c*^-order Euler-Poincare reduction). Let L : T^'^^G be a G -invariant La- 

grangian and let i : kg ^ K. be the associated reduced Lagrangian. Let g(t) be a curve in G and 
^{t) = g{t)g{t)~^, resp. ^(t) = g{t)^^g(t) be the reduced curve in the Lie algebra q. Then the 
following assertions are equivalent. 

(i) The curve g{t) is a solution of the k^^-order Euler-Lagrange equations for L : T^^^G — ?► M. 

(ii) Hamilton's variational principle 

5 L{g,g,...,g'^''^)dt = 

holds upon using variations 5g such that 5g^^^ vanish at the endpoints for j = 0, k — 1. 

(iii) The k^^-order Euler-Poincare equations for i : kg M.- 

fc— 1 

{dt±s.dl)J2i-iyd{^ = 0. (3.6) 

j=0 ^ 



(iv) The constrained variational principle 



holds for constrained variations of the form 6C, = dtrj =]= [C,,'']], where rj is an arbitrary curve 
in Q such that r]^^^ vanish at the endpoints, for all j = 0, ...,k — 1. 
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Remark 3.2. The quotient map (2.14), respectively (2.15), can be used for any Lie group G. In 
the case of matrix groups, one might consider the alternative quotient map of the form 

{9,9,-;9^^^) ^ ii^i,--,i^k) , i^j := g^^'^ respectively uj := g'^g^^\ (3.7) 
One may easily pass from the variables ^, dl^~^^^^ to the variables (i^i, Uk)- For example: 

i = dt{gg~^) = gg~^ - gg~^gg~^ = z/2 - i^ii^i (3.8) 

^ = Z/3 - 2z/2Z/l + 2z/iZ/iZ/i - Z/iI/2, 

and so forth, by using the rule uj = Uj^i — UjUi. Here all concatenations mean matrix multiplica- 
tions. One can easily derive the constrained variations and the /c*'*-order Euler-Poincare equations 
associated to this quotient map in a similar way as above. 



3.2 Example: Riemannian cubics 



In this section we apply the k^^-oidei Euler-Poincare reduction to the particular case of 2-splines 
on Lie groups. Fix a right-, respectively left-invariant Riemannian metric 7 on the Lie group G. 
We denote by 

the corresponding squared norm of a vector Vg G TgG. The inner product induced on the Lie 
algebra g is also denoted by 7 : g x g — )• M and its squared norm by 



We recall that the associated isomorphisms 

b : ^ 0*, ^^t, and 

are defined by 



0^0, /i /i^ 



{^\v) =li^,v), for all e,r7G0, and ^ := \) \ 
where ( , ) denotes the dual pairing between q* and q. 

Proposition 3.3. Consider the Lagrangian L : T^'^^G — )■ M for geometric 2-splines, given by 

2 



(3.9) 
(3.10) 



L{g,g,g) 



D 
Dt 



(3.11) 



where \\ ■ \\ is the norm of a right-, respectively left-invariant metric on G. Then L is right-, 
respectively left-invariant and induces the reduced Lagrangian £ : 20 — )■ M given by 



i{^,0 = l e±adje"' 



(3.12) 



where ad^ by adlrj := (ad^(r7''))^, for any ^,77 G 0. 
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Proof. Let us recall the expression of the Levi-Civita covariant derivative associated to a right 
(respectively left) G-invariant Riemannian metric on G. For X G X{G) and Vg e TgG, we have 
(e.g., [KM97], Section 46.5) 

V.Xig) = TRg (^dfivg) + ^ adt fig) + i ad^^^^ v-^[v, fig)]^ , v := v^g-' (3.13) 



(111 
^fivg) - 2 fig) - - ad|(^) + - [ 



vJig)] 



v:= g Vg 



(3.14) 



where / G ^-"((7; g) is uniquely determined by the condition X{g) = TRg{f{g)) for right-, respec- 
tively X{g) = TLg{f{g)) for left G-invariance. Therefore, we have 



D 



git) =Wgg = TRgU + -^d\i+- ad^^ - -[i.i] = TRg U + ad^^ 



respectively 



D 
Dt 



git) = Vgg = TLg (^^ - 1 ad^ ^ - ^ ad^ ^ + ^ [^, = TLg - adj , 



99' 



where we used Xig) = g, Vg = g, so fig) 
dfivg) = t 

Thus we obtain, due to the right-, or left-invariance of the metric 7 

2 



^ (respectively, fig) = g ^g = ^) and 



Li9,9,9) 



1 



D 
Dt 



e±adje 



(3.15) 



which depends only on the right invariant quantity ^ = gg ^, respectively the left invariant quantity 
^ = g^^g. Accordingly, L is right-, or left-invariant, and the group-reduced Lagrangian is 



e±adU 



□ 



which completes the proof. 

Remark 3.4. The above considerations generalize to geometric /c-splines for k > 2. Indeed 
iterated application of formulas (3.13), (3.14) yields 

j^g = TRg ir]k) , respectively —g = TLg (r^^) , 
where the quantities rji^ ^ g are defined by the recursive formulae 



r]i=i± adl ^, and r]k = rjk-i ± ^ (adj r]k^i + ad],^,_^ ^ + ad^,_j 



(3.16) 



for ^ = gg ^, respectively = g ^g- Therefore, the Lagrangian (2.9) for geometric A;-splines on a 
Lie group G with right-, respectively left-invariant Riemannian metric, 

2 



Lk {g,g,...,g^^>) 



1 



D 



k-l 



Dt^- 



-i9 



is right-, respectively left-invariant, and the reduced Lagrangian is 



^(e,e,...,e 



(fc-l)N 



2 11%-lllg 



(3.17) 
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Computing the second-order Euler-Poincare equations for splines. Let us compute the 
Euler-Poincare equations for k = 2. The required variational derivatives of the reduced Lagrangian 
(3.12) are given by 

^^^=e±^die=--v' and | = ^(ad;e^ + (ad,0') efl*. (3.18) 

2"'^-order Euler-Poincare equation 
= 0, with 7]^:=t±adlt, (3.19) 

= 0, with r/:=^±adJC. (3.20) 



These are the reduced equations for geometric 2-splines associated to a left-, or right-invariant 
Riemannian metric on the Lie group G. 

In an analogous fashion one can derive the Euler-Poincare equations for geometric fc-splines, 
using the reduced Lagrangian (3.17). 

When the metric is left-, and right-invariant (bi-invariant) further simplifications arise. 



From formula (3.5) with k = 2 one then finds the 

(9t±ad*) (dtv'±^d;t±{^d,0') 
or, equivalently, 

feiadj) (at77±adi,e±ad^e) 



Example 1: Bi-invariant metric and the NHP equation. In the case of a bi-invariant 
Riemannian metric, we have ad^ f] = — ad^ rj and therefore i]^ = C^ so that f] = ^ and the equations 
(3.20) become 

{dt ± ad*) = or (dt ± adj) ^ = or T ^] = 0, (3.21) 

as in [(^895]. Note that in this case, the reduced Lagrangian (3.12) is simply given by = 
^H'CP- We also remark that since the metric is bi-invariant, one may choose to reduce the system 
either on the right or on the left. This choice will determine which sign appears in (3.21). 

Taking G = SO {3), we recover the NHP equation (1.1) of [NHP89]: 

!r2 = ±rixfi. (3.22) 

In [NHP89], the unreduced equations in the general case are also derived, but the symmetry 
reduced equation is given only for SO (3) with bi-invariant metric. 

Remark 3.5. [Conventions for so(3) and so(3)*] 

In equation (3.22) and throughout the paper we use vector notation for the Lie algebra so(3) of 
the Lie group of rotations 50(3), as well as for its dual 50 (3)*. One identifies so (3) with M.^ via 
the familiar isomorphism 

/ a\ / -a b \ 

^:M3->so(3), n=lb\^n=n=l a O -c , (3.23) 

V c ; \-b c J 

called the hat map. This is a Lie algebra isomorphism when the vector cross product x is used as 
the Lie bracket operation on M^. The identification of so (3) with induces an isomorphism of 
the dual spaces so(3)* = (M^)* = 
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Example 2: Elastica. Another example of the 2"'^-order Euler-Poincare equation arises in the 
case of elastica treated in [HB04b], whose Lagrangian is 



T 



mi; + 2 



D 

m 



and whose reduced Lagrangian is 
Using the 2"'^-order Euler-Poincare equation (3.5) one easily obtains the reduced equations 



(3.24) 



(dt ± ad|) {dtri ± adj ^ ± ad^ ^ - r^^) = 0, with r/ := ^ ± ad^ ^, 



(3.25) 



which simplify to 



in the bi-invariant case. 



[d, ± adt) {d^^ - r^O = 



Remark 3.6. We now consider the particular case G = 5*0(3). Let I be a 3 x 3 symmetric positive 
definite matrix (inertia tensor) and consider the inner product 7(rii, 0.2) = If^i ■ ^2 on R^. The 
Lagrangian for the elastica on 5*0(3) reads 



L(A,A,A) 



where || ■ ||a is the right-, respectively left-invariant metric induced the inner product 7. Relative 
to this inner product we have 

so the reduced Lagrangian (3.24) reads 



r2 


A 


2 1 










y 




» + 2 





~2 



fi-ifi + ^ (in ±inxny (in ±mxny 



(3.26) 



If r = 0, this expression can be interpreted as the Lagrangian for geometric 2-splines of a rigid 
body. 

If I is the identity, the Lagrangian in (3. 26) simplifies to 

1 . . 

^ ' ' 2 2 



and the Lagrangian of the NHP equation is recovered when r = 0. 
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Example 3: L^-splines. One can consider geometric 2-splines on the diffeomorphism group 
of a manifold V as follows. Fix a Riemannian metric g on V and consider the associated 
right-invariant Riemannian metric on G = Diff(D) and its induced second-order Lagrangian 

Hv,v,v) = \ 

on T*^^^ Diff (P). The reduced Lagrangian on 20 = 2X('D) reads 

^(u,u) = ^||u + adjjuf , 

where ad^ denotes the transpose with respect to the inner product, given by adj,v = VuV + 
(Vu)""" ■ V + vdivu. In this case, the ad^ and ad terms in (3.20) combine to produce the spline 
equation 

ipt + ad^) {dtv + 2Sv ■ u + u div v) = 0, v = S^u + 2Su ■ u + u div u, 

where Su := (Vu + Vu"^) /2 is the strain-rate tensor. 

In the incompressible case, that is when G = Diff^oi(T'), the transpose of ad relative to the 
inner product on divergence free vector field is denoted by ad"*" and is related to ad^ by the formula 

ad+ V = P (adj, v) = P (^VuV + (Vu)"^ ■ v) , 

where P denotes the Hodge projector onto the divergence free vector fields. In this case (3.20) 
reads 

{dt + ad+) (^t^r + P (v^u + (Vv)^ ■ + VuV - Vvu) =0, v = dtU+2F (Su ■ u) , div u = 0. 

Remarkably, using the formula adjj Vp = V(Vp ■ u) for divu = 0, all the gradient terms arising 
from the Hodge projector can be assembled in a single gradient term, thereby producing the 
incompressible 2-spline equations 

[dt + adjj) {dtw + 2Sw ■ u) = — Vp, w = dtu + 2Su ■ u, div u = 0, 

where ad^ (and not ad"*") is used. 



D 

Dt" 



Example 4: i/^-splines. One can alternatively consider splines relative to the right-invariant 
metric induced by the inner product {Qu, u), where Q = {1 — a^A). In this case, the 2-spline 
equation reads 

{dt + ady {dtQ^r + 2S{Q ■ adv) • u) = 0, Qv = OtQu + ad^, Qu (3.27) 
where S(L) := 1{L + L*). 

Remark 3.7. Note that in order to obtain the simple expression 
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(instead of (3.24)) for i by reduction, one needs to modify the spline Lagrangian as 

L{9,9,9) = ^11^11? + ^ 

where || • i = 1,2, are two norms associated to two G-invariant Riemannian metrics on G and 
ad^ is extended as a bihnear map TgG x TgG — t- TgG by G-invariance. 

The associated 2"'^ order Euler-Poincare equations are simpler than the one associated to 2- 
splines. For example, the reduced Lagrangian £(u, li) = ^ ((1 — a^A)u, u) + ^Hup produces the 
following modification of the EPDiff equation: 

{dt + ady (u - (a^ Au + Mu)) = 0. (3.28) 
3.3 Parameter dependent Lagrangians 

In many situations, such as the heavy top or the compressible fluid, the Lagrangian of the system 
is defined on the tangent bundle TG of the configuration Lie group G, but it is not G-invariant. In 
these cases, the Lagrangian depends parametrically on a quantity go in ^ manifold Q on which G 
acts and that breaks the symmetry of the Lagrangian L = Lg^. We refer to [HMI198], [GBR09] for 
the case of (affine) representation on vector spaces, relation with semidirect products and many 
examples. This theory was extended to arbitrary actions on manifolds in [GBTIO] for applications 
to symmetry breaking phenomena. We now briefiy present the extension of this theory to the case 
of higher order Lagrangians. 

Consider a fc*''-order Lagrangian : T^^^G R, depending on a parameter go in a manifold 
Q. We suppose that G acts on the manifold Q and that the Lagrangian L is G-invariant under 
the action of G on both T^'^'^G and Q, where we now see L as a function defined on T^'^^G x Q. 
Concerning the action of G on T^'^^G x Q, there are several variants that one needs to consider, 
since they all appear in applications. 

(1) First, one has the right, respectively the left action 

{g,g,...,g^^\qQ) ^ {gh, gh, g^''^h,h-\o) , 
respectively {g, g, g^''\ qo) {kg, kg, hg^^\ qoh-^) . 

The reduced variables are (.^,g) = {gg~^, gqo), respectively (^, g) = {g^^ 9 , log) ■ In this case, the 
A;*'*-order Euler-Lagrange equations for L^^ on T^'^^G (where go G Q is a fixed parameter) are 
equivalent to the fc^'^-order Euler-Poincare equations together with the advection equation 



{9t ± ad*) \T,{-iydi^] = ^ (^) ' ~ ^'^^^^ = °' ^^-^^^ 



with initial condition go. Here ^^(g) denotes the infinitesimal generator of the G action on Q and 
J : T*Q — i- g* defined by (J(aq),^) := (ag,^Q(g)) denotes the momentum map associated to the 
G action on T*Q. 



— q ± adt q 
Dt^ 9^ 
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The associated variational principle reads 

fi2 



5 J\{^,i,...,^^'-'\q)dt = 0, 



relative to the constrained variations (3.3) and constrained variations of q given by 5q = riQ{q), 
where 77 = {Sg)g~^, respectively rj = g~^[6g). Equations (3.29) and their variational formulation 
can be obtained by an easy generalization of the approach used in Section §3.1. 

(2) Secondly, one can consider the right, respectively the left action 

{g,g,...,g^''\qo) ^ {gh, gh, g^^^h, qoh) , 
respectively {g, g, g^''^ qo) ^ {hg,hg, ...,hg'^''\hqo) . 

The reduced variables are {C,,q) = (fi'5'~^, Q'o5'~^)) respectively {C,,q) = {g^^g, g^^Qo) and one gets 
the reduced equations 



{d, ± ad*) (jZ^-iy^i^^ = -J (^) , d,q + ^Q{q) = 0, (3.30) 



with initial condition qo. 

U Q = V* is the dual of a G-representation space so which G acts on V* by the dual represen- 
tation, the above equations reduce to 

{dt±-di)(j2yydi^^=f^<>a, (3.31) 

where the diamond operation o : V x V* ^ Q* is defined by 

{voa,0 = {a,^v{v)) , 

and therefore J(a, v) = —voa. These are the higher order version of the Euler-Poincare equations 
with advected quantities studied in [HMR98]. 

Example : Rate-dependent fluid models. Rate-dependent fluid models are usually defined 
using Lagrangians that depend on the strain-rate tensor S := (Vu + (Vu)"'")/2 and its higher spatial 
derivatives [BFHL88]. A related class of spatially- regularized fluid models have been introduced 
as turbulence models [FHTOl]. 

Yet another class of rate-dependent fluid models may be defined, e.g., as 2nd order Euler- 
Lagrange equations T*^^^ DifF(T') for a parameter dependent Lagrangian La^- The group reduced 
representation of the equations of motion for such rate-dependent fluids is found from the previous 
manipulations, namely 

(dt±£u)(^-dt—] =^oa, idt + £u)a = 0. (3.32) 
\ ou ou / da 
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One has the + sign for right and the — sign for left invariance. The usual Eulerian fluid represen- 
tation is right-invariant and so takes the + sign. Physically, these fluid models penalize the flow for 
producing higher temporal frequencies. Therefore, these models might be considered as candidates 
for frequency-regularized models for fluid turbulence. The Kelvin theorem for these fluids involves 
circulation of the higher time derivatives. For right-invariant higher-order Lagrangians, the Kelvin 
theorem becomes 

d r I f 5i ^6i\ r 1 6i 



dt /c{u) L) \5u Su J /^(u) D 5a 

where the density D satisfies the continuity equation {dt + £u)D = 0. Consequently, the integrands 
in the previous formula are 1-forms and thus may be integrated around the closed curve c(u) 
moving with the fluid velocity, u. This is the statement of the Kelvin-Noether theorem [HMR98] 
for fc-splines. 

3.4 Splines with constraints 



Suppose that one wants to minimize the action J^^ L{q, q, q)dt over curves q{t) E Q subject to the 
condition 

{u^{q),q) = ki, i = l,...,fc, 
where Ui are 1-forms and ki G M. One uses the variational principle 

S ^ ' (^L{q, ^' ^) + 5^ ^) - ^^)) = 0' (3.33) 

for arbitrary variations 6Xi of the curves Xi(t), i = 1, k, and for variations 6q vanishing at the 
endpoints. Variations relative to q yield the equation 

^2 dL ddL dL ^ 



i=l 



dt"^ dq dt dq dq 

whereas variations relative to Aj yield the constraint. For example, for the Lagrangian (2.8) this 
yields the equations 



m + R I —q{t). m ) m = -^^^M) + (Aa^du;. + a, 
^ ^ i=i 



as in [BC93], [CS95], [HB04a], where Ui = X\ for given linearly independent vector fields Xi, . . . G 
X(Q). 

Remark 3.8. The variational principle (3.33) is equivalent to 

5 ^L(g,g, g) + ^ Ai (a;i(g),g)^ (it = and {uji{q),q) = ki, i = l,...,k, (3.34) 
where only variations of q are involved and the term containing k^ is suppressed. 
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We now consider the special case Q = G and we suppose that the one-forms cjj, i = 1, k, 
are G-invariant. That is, we can write 

{uJi{g),Vg) = {uji{g)g-'^, ^gg'^) = {Q, Vgg~^) resp. {uJi{g),Vg) = {Q, g'^Vg) , 

where 

:= uji{e) e 0*. 

The reduction of the variational principle (3.33) yields the constrained variational principle 



for arbitrary variations 5\i of the curves Aj(t), i = l,...,k, and variations of ^{t) satisfying the 
constraints (3.3). Equivalently, using (3.34) we rewrite the stationarity condition as 



6 J^' (£{^,O + {z,o)dt = and {C^,0 = k^, ^ = 1,...,A;, 



for variations of C, satisfying (3.3) and where we have defined z := Yli=i ^id ^ 0*- We obtain the 
equations 



{dt ± ad^) i^- - 9,^ + ^ j = 0, (0, = h. 



With the Lagrangian (3.12) for 2-splines, we find the reduced equations 

{^t ± ad^ j [dtf] ± ad|^ ^ ± ad^ C ~ -2) =0, with 77 := ^ ± ad^ ^ 

and for bi-invariant metrics, we get 

(a^iadj) (e-z) =0, i.e., e T[e,e-^]-i = 0, {Ci.i) = h, 

which coincides with equation (39) in [CS95]. See also [BC96a]. 

We can also consider higher-order constraints, with associated variational principle 



5 



to 



^{i.L ...,^'^'-'^) + (^0,0 + - + {zk-i,^^'-''>)) dt = 0. 



In this case, one obtains the equations 

fc-i 



(9.±adax:(-i)^a^(^+-.)=o. 



i=o 

For example, with k = 2, and a bi-invariant metric we have 



(dt ± adj) + zo) = 0. 
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4 Clebsch-Pontryagin optimal control 

Here we develop the /c*''-order Euler-Poincare equations from an optimal control approach. The 
ideas in [GBR 10] for k = 1 generalize easily to higher order. 

Definition 4.1. Let ^ be a (right or left) action of a Lie group G on a manifold Q. For a Lie 
algebra element ^ E Q let 

'■= 37 ^exp{to(5)' 

denote the corresponding infinitesimal generator of the action. Given a cost function £ : /eg — )■ M, 
the Clebsch-Pontryagin optimal control problem is, by definition, 



mm 

m 



£ i[^,i....,^^''-'^)dt (4.1) 



subject to the following conditions: 

(A) q = Uq) or (A)' q = -^Qiq); 

(B) g(0) = go and q{T) = qx; 

(C) i^^\Q)=ii and e'KT)=eT, J=0,...,k-2, 
where qo, qx & Q and ^q, G 0, j = 0, /c — 2, are given. 

Variational equations. We suppose that condition (A) of Definition 4.1 holds. (The calculation 
for case (A)' is similar.) The resolution of this problem uses the Pontryagin maximum principle 
which, under sufficient smoothness conditions, implies that its solution necessarily satisfies the 
variational principle 

S [i[^,i-;^^'-'^) + {<y,q-^Qiq)))dt = 0, 

for curves t h-> ^(t) G 5 and t 1— )■ a{t) G T*^^-^Q. This variational principle yields the conditions 

fc— 1 

J(«W) = E(-l)'^*iu) a = iT*Q{al (4.2) 

in which J : T*Q — )■ g* is the cotangent bundle momentum map, as above, and ^t*q denotes the 
infinitesimal generator of the cotangent lifted action, denoted : G x T*Q — )■ T*Q. 

If G acts on the right (respectively left), a solution a{t) of d = ^t*q{ol) is necessarily of the 
form a{t) = <l>Jj;)(a(0)), where g{0) = e and g{t)g{t)-^ = ^(t) (respectively g{t)-^g{t) = ^(t)). 
The above conditions imply coadjoint motion, 

3{a{t)) = J ($ J(;)(«(0))) = Ad;(,) J(«(0)), respectively J(a(t)) = Ad;^,^, J(«(0)), 

and by differentiating relative to time, we obtain the left (right) Euler-Poincare equations: 

^J(a(t)) = ad*(4)-i^(t) J(a(t)) = ad|(^) J(a(t)), 
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respectively 

— J(a(t)) = -ad*(i)3(i)-i J(a(t)) = - ad^(i) J(a(t)). 
Upon using the first condition in (4.2), we recover the k^^-oidei Euler-Poincare equations 

fc-i 



idtT^dl)J2i-^yd{^ = 0. (4.3) 

j=0 



Example 1: Clebsch approach to the NHP equations. The NHP equations can be obtained 
from the Clebsch approach by considering the action of 50(3) on M^. The Clebsch- Pontryagin 
control problem is 



T 

2, 



min / ||ri|pdt, subject to q = i7 x q, q(0) = qo, q(T) = qr, ^1(0) = ^o, ^(T) = ^t- 







The stationarity conditions (4.2) read q x p = — f2, q = f2 x q and p = $7 x p. One directly 
observes that they imply the NHP equations. 

Example 2: Clebsch approach to if ^-splines. We let the diffeomorphism group Difft,o;(T') 
act on the left on the space of embeddings Emb(S', "D) of a manifold S in D. The associated 
Clebsch-Pontryagin control problem is 



mm 



in / ||u + adjj u||^i (it, subject to Q = u o Q, Q(0) = Qo, Q(T) = Qt, 
*) Jo 

u(0) = uo, u(T) = vlt- 



The condition 

J(Q, P) = dtQ^ + 2S(g ■ adv) ■ u, 

where Qv = dtQ\i + adJjQu, together with the Hamilton equations on T*Emb(5', P) imply the 
spline equations (3.27). In the case of equations (3.28) the condition (4.2) reads 

J(Q,P) =Qu-Ui,. 



Additional g-dependence in the Lagrangian. One can easily include a g-dependence in 
the cost function £ of the Clebsch-Pontryagin optimal control problem (4.1). In this case, the 
stationarity conditions (4.2) become 



k-l 



J("(^)) = 5^*^"^)^^*^ d = ^r-Q(a)+Ver„— , 



(4.4) 



where for a, /3 G T*Q, the vertical lift of /3 relative to a is defined by 



Ver„/3 : = 



d 
ds 



{a + sf3) eT^{T*Q). 



s=0 
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In this case, the differential equation (4.3) for the control ^{t) generalizes to 

(a,^adj)|:(-l)'a?-^=j(|), (4.5) 

where 3 : T*Q ^ g* is again the cotangent bundle momentum map. 
Remark 4.2. [Recovering Euler-Poincare equations of §3.3] 

Equations (4.5) for the control recover the fc^^^-order Euler-Poincare equations (3.29). Note that 
a right, respectively left action of G on Q produces the left, respectively right Euler-Poincare 
equations in (4.3) consistently with the results in §3.3. In order to obtain the Euler-Poincare 
equations (3.30) one needs to impose condition (A)' instead of (A) on the dynamics on q. 

5 Higher-order template matching problems 

In this section we generalize the methods of [BGBHRIO] to higher order because the added smooth- 
ness provided by higher-order models makes them attractive for longitudinal data interpolation, 
in particular in Computational Anatomy ( CA ). 

We first give a brief account of the previous work done on longitudinal data interpolation in CA. 
Then we derive the equations that generalize [BGBHRIO]. After making a few remarks concerning 
the gain in smoothness, we provide a qualitative discussion of two Lagrangians of interest for CA. 
Finally, we close the section by demonstrating the spline approach to template matching for the 
finite dimensional case of fitting a spline through a sequence of orientations on SO {3). 

5.1 Previous work on longitudinal data interpolation in CA 

CA is concerned with modeling and quantifying diffeomorphic evolutions of shapes, as presented 
in [MTY02, MYOl]. Usually one aims at finding a geodesic path, on the space of shapes, between 
given initial and final data. This approach can be adapted for longitudinal data interpolation; that 
is, interpolation through a sequence of data points. One may interpolate between the given data 
points in such a way that the path is piecewise-geodesic, [BK08, DPG+09]. It was, however, argued 
in [TVIO] that higher order models, i.e., models that provide more smoothness than the piecewise- 
geodesic one, are better suited as growth models for typical biological evolutions. As an example 
of such a higher-order model, spline interpolation on the Riemannian manifold of landmarks was 
studied there. In the next paragraph, we will consider another class of models of interest for CA 
that are inspired by an optimal control viewpoint. Indeed, the time-dependent vector field can 
be seen as a control variable acting on the template and the penalization on this control variable 
will be directly defined on the Lie algebra. Finally, we underline that this class of models is an 
interesting alternative to the shape splines model presented in [fVlO]. 

5.2 Euler-Lagrange equations for higher-order template matching 

Let G be a Lie group with Lie algebra q, and let 



GxV^V, {g,l)^gl 



(5.1) 
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be a left representation of G on V. Let || • be a norm on V. We consider minimization problems 
of the following abstract form: 

Given a Lagrangian £ : {k — l)g R, cr, ti, . . . , t/ G M, Tq, Jj^, . . . , Jj, e V, and^^...,^^ 2 g 
minimize the functional 

E[^] ■■= r mt), • • • , ^'•'-'\t))dt + ^ E - (5-2) 

i=i 

subject to the conditions .^'•■'•'(O) = ^q, j = 0, . . . , k — 2, where g^{ti) is the flow of ^(t) evaluated at 
time ti. The minimization is carried out over the space 

Vu-i := G C'-\%t,U) I ^'^'-'^ e C^^^ (([t,„tm]):=i)} 

where (([tj, tj+i])'^Q) denotes the set of piecewise smooth curves whose only discontinuities 
would be at the ti, i = 1, . . . ,1 — 1, i.e. 

C^c^ (([t.,tm]):=S) ■= {/ e ^'([0,T],s) I = 0, . . . , / - 1 3r G C-([t„ t,+i], g) s.t. / = /|^(,^,^^^)] 



More precisely, given such a curve ^(t) in the Lie algebra g, its flow g^ : t ^ 9^{t) G G is a 
continuous curve defined by the conditions 

g^{0) = e, and |/(t) = , (5.3) 

whenever t is in one of the open intervals (0, ti), . . . , t/). Here we used the notation ^(t) (7^ (t) := 
TRg^t)^{t). 

We typically think of [It-,^, . . . , It^) as the time-sequence of data, indexed by time points tj, 
j = 1,.../, and To is the template (the source image). Moreover, ^ : t G [0,^;] i-> C,{t) G 
is typically a time-dependent vector field (sufficiently smooth in time) that generates a flow of 
diffeomorphisms g^ : t & [0,^/] ^ g^{t) ^ G. Note that, in this case, the Lie group G is infinite 
dimensional and a rigorous framework to work in is the large deformations by diffeomorphisms 
setting thoroughly explained in [TY05]. We will informally refer to this case as the diffeomorphisms 
case or infinite dimensional case. The expression g^{ti)TQ represents the template at time tj, as it 
is being deformed by the fiow of diffeomorphisms. Inspired by the second-order model presented in 
[TVIO], this subsection thus generalizes the work of [BGBHRIO] in two directions. First, we allow 
for a higher-order penalization on the time-dependent vector field given by the first term of the 
functional (5.2); second, the similarity measure (second term in (5.2)) takes into account several 
time points in order to compare the deformed template with the time-sequence target. 

Staying at a general level, we will take the geometric viewpoint of [BGBHRIO] in order to 
derive the Euler-Lagrange equations satisfied by any minimizer of E. We suppose that the norm 
on V is induced by an inner product {,)y and denote by b the isomorphism 

that satisfies 

(/, J)y = (/^ J) for all I,J eV, 
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where we wrote ( , ) for the duahty pairing between V and its dual V*. The action (5.1) of G on 
V induces an action on V*, 

GxV*^V*, {g,u)^gu 

that is defined by the identity 

{goo, I) = (w, g-^l) for all I e V, 00 e V* , g e G. (5.4) 

The cotangent-hft momentum map o : V x V* ^ Q* for the action of G on V is defined by the 
identity 

{lou,0 = {uJ,^I), for all I eV, u eV*, ^ eg, (5.5) 



where the brackets on both sides represent the duality pairings of the respective spaces for q and 

— I 

dt 1 1=0 ■ 



V, and where denotes the infinitesimal action of q on V, defined as ^ / := 4:1 „g{t)I G V for 



any G^ curve g : [— — )■ G that satisfies g{0) = e and ^l^^^^git) = ^ G 0. Note that equations 
(5.4) and (5.5) imply 

Adl-i{Iouj) = glogco. (5.6) 
For the flow defined in (5.3), we also introduce the notation 

gl ■■= gHt) {9Hs)y' . (5.7) 

Lemma 2.4 in [BGBHRIO], which is an adaptation from [Via09] and [BMTY05], gives the 
derivative of the flow at a given time with respect to a variation {e,t) ^ ^{t) + eSC,{t) G of a 
smooth curve ^ = ^o- Namely, 



de 



9t = 9is f (Ad^|/e(r))rfrGT , G. 

=0 J s ^ ' 



(5. 



Importantly, formula (5.8) also holds for the diffeomorphisms case in a non-smooth setting, as 
shown in [TY05], where the assumption is ^ G L^([0, t;], \^). Moreover, this proof can be adapted 
to the case of a finite dimensional Lie group. In particular, formula (5.8) can be used for ^ G Vk~i, 
whether one works with finite dimensional Lie groups or diffeomorphism groups. 

Formula (5.8) and equation (5.6) are the key ingredients needed in order to take variations of 
the similarity measure in (5.2). With these preparations it is now straightforward to adapt the 
calculations done in the proof of Theorem 2.5 of [BGBHRIO] to our case, in order to show that 
the following theorem holds. 

Theorem 5.1. A curve C, G Vk-i is an extremal for the functional E, i.e., 5E = if and only if 
(I), (II), and (III) below hold: 



(I) For t in any of the open intervals (0, ti), . . . , (t^-i, ti), 

k—l I 

E(-l)^'^5§) =-J2^i0Mit) (gioToogly) , (5.9) 
where vr* is defined by 



and Xlo,u] is the characteristic function of the interval [0,ti]. 



TT : = 



Gay-Balmaz et al. 



Lie group reduction of higher-order invariant variational problems 



26 



(II) For z = 1, 1 andr = 0,...,k-2, 



lim y (-ir^-^^^ — —At) = lim y {-iy-^-^^—^{t). (s.io) 

' j>r+l ^ » i>'r+l 



(III) Forr = 0,...,k-2, 



j>r+l ^ 



Note that there is no condition at to = analogous to (III) because of the fixed end point 
conditions ^(^^(0) = for j = 0, . . . , A; - 2. 

Proof. Set to = for convenience. A series of partial integrations taking into account the fixed 
end point conditions ^'•■''^(0) = $,1, j = 0, . . . , k — 2, leads to 



i=0 \ j=0 

I— I k—2 I I fc— 1 / To—r — 1 Ai—'r—\ 



i=l r=0 \ \i>r+l ^ ^ ^ y 

+(E(-ir'-'|S^(«<).^«"(«<)))- (5.12) 

Note that the hypothesis ^ G Pfc-i is sufficient to give meaning to the previous formula. 

On the other hand, using formula (5.8) and mimicking the computations done in [BGBHRIO], 
one finds for the variation of the similarity measure 

^ (^^E - = \l (i2xio,tAt) {gioToogly),6at)^ dt . (5.13) 

Assembling the two contributions to 6E, we arrive at 

6E = y^ J|) (^) + E^[o.*.]W {9ioTooglx),6m) dt 



+EE f( E (-ir-' f^^(*.-) - S^^jyC?) ) 



i=l r=0 \ \i>»'+l 



+E(E(-iy-'-'|;i^(«,^?'"(«)). (5.14) 

r=0 \i>r+l ^ / / 

Stationarity bE = therefore leads to equations (5.9) - (5.11). □ 
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Remark 5.2. The right-hand side of equation (5.9) follows coadjoint motion on every open interval 
(0, ti), . . . , {ti, t;_i). Therefore, 

(| + «dj,„)|:(-l)^-|^^0. (5,15) 
in which we once again recognize the higher-order Euler-Poincare equation (3.4). 



5.3 Two examples of interest for computational anatomy 

Regarding potential applications in CA, an interesting property of higher-order models is the 
gain in smoothness of the optimal path T : t E [0,ti] t— )■ g^{t)To e V, in comparison with first- 
order models. For instance, in the case of piecewise-geodesic (i.e., first-order) interpolation, where 
■= IW^WI, equation (5.9) reads 

I 

^(t) = -Y.^lo,u]it) [gioToOglx) . (5.16) 

i=l 

In general therefore, ^ will be discontinuous at each time point tj for i < I, which implies non- 
differentiability of T at these points. In contrast, for the Lagrangian := ^H^Hg, equation (5.9) 
becomes 

m = J2^ioMit) {gioTooglx) ■ (5.17) 

Now the curves ^{t) and T{t) are and on [0,t;], respectively. Note that the inexact inter- 
polation we consider here yields a curve T, whereas the exact interpolation method presented 
in the example of Section 2.2 leads to solutions. Note also that the minimization of the 
functional E for ii when / = 1 produces Lie-exponential solutions on G. More precisely, if the 
Lie-exponential map is surjective and the action of G on is transitive, then there exists ^ 
such that exp(ti^o)^o = ^ti- Hence, the constant curve ^ = is a minimizer of the functional E, 
with E[C] = 0. The Lie-exponential has been widely used in CA, for instance in [AFPA06, Ash07]. 

Another Lagrangian of interest for CA is £2(C, := lU + adj ^llg' ^hich measures the acceler- 
ation on the Lie group for the right-invariant metric induced by the norm || ■ ||g. The Lagrangian 
£2 may therefore have more geometrical meaning than ii. However, ii is worth studying since 
it is simpler from both the computational and the analytical point of view: The existence of a 
minimizer for ii can be obtained straightforwardly following the strategy of [TY05]. In contrast, 
a deeper analytical study is required for £2, since analytical issues arise in infinite dimensions. 



5.4 Template matching on the sphere 

Consider as a finite-dimensional example G = S0{3) with norm ||f2||go(3) = y/Cl ■ Ifi on the Lie 
algebra so (3), where I is a symmetric positive-definite matrix (the moment of inertia tensor). Let 
^ = with II ■ ||]fj3 the Euclidean distance. We would like to interpolate a time sequence of points 

on the unit sphere S'^ C M'^ starting from the template Tq = . Choose the times to be 
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ti = |i for z = 1, . . . , 5, and 

= I 1 I , = I I , = — I I , J,^ = — I 1 I , I,^ ^ 





(5.18) 



The associated minimization problem for a given Lagrangian . . . , fl^^ is: 
Minimize 

E[n] := em), n('^-'\t))dt + ^ E ll^"(^Oro - lu ||k3 , (5.19) 
subject to the conditions fl^^\0) = fig, j = 0, . . . ,k — 2, where A"(t) is a continuous curve defined 

by 

A"(0) = e , and ^A"(t) = n{t)A''{t) , 

whenever t is in one of the open intervals (0,ti), . . . , (^4,^5). As we mentioned in Section 5.3, an 
important property of higher-order models is the increase in smoothness of the optimal path when 
compared with first-order models. We illustrate this behavior in Figures 5.1 and 5.2: 

Figure 5.1 shows the interpolation between the given points J^^, . . . , J^g for the first order La- 
grangian 

= ■ in. (5.20) 

We contrast this with the second order model 

n) = ^(n + r\n x m)) + r\n x m)) . (5.21) 

Note that this is the reduced Lagrangian for splines on 5*0(3), as we discussed in Section 3, and 
for I = e we recognize equation (5.15) to be the NHP equation (3.22). 

Figure 5.2 visualizes the resulting interpolation for two different choices of the moment of inertia 
tensor I, namely 

/100\ 1 f ^ ^ ^\ 

Ii := 1 and I2 := ^ 2 . (5.22) 

In order to compare the two cases we have normalized I2 in such a way that it has the same 
norm as Ii with respect to the norm ||I|p = tr(l"'"l). The figures were obtained by minimizing 
the functional E using the downhill simplex algorithm fmin_tnc that is included in the optimize 
package of SciPy, [.T()P+ ]. 
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(a) a = 0.18 (b) a = 0.01 



Fig. 5.1: First order template matching results are shown for the Lagrangian (5.20) with I = e, for two different 
values of tolerance a. These values have been chosen so that the sum of the mismatch penalties is similar in size 
to the one obtained in the second order template matching shown in Figure 5.2. As might be expected, when the 
tolerance is smaller, the first order curves pass nearer their intended target points. These first order curves possess 
jumps in tangent directions at the beginning of each new time interval. 




(c) a = 0.05, I = h (d) (T = 0.001, I = h 



Fig. 5.2: The pictures in the top row show the template matching for the Lagrangian (5.20) with Ii with two 
different values of tolerance, a. The bottom row represents the corresponding matching results for I2. One observes 
that the quality of matching increases as the tolerance decreases. This is due to the increased weight on the penalty 
term in (5.2). The color of the curves represents the magnitude of the velocity vector of the curve on the sphere 
(red is large, white is small). We fixed the initial angular velocity f2(0) = ^ (0, 0, 1). On comparing these figures 
with those in the first order case, one observes that the second-order method produces smoother curves. 



Gay-Balmaz et al. 



Lie group reduction of higher-order invariant variational problems 



30 



Remark 5.3. Standard variational calculus arguments ensure the existence of a minimizer to 
to the functional (5.2) with Lagrangian (5.21). In Theorem 5.1, we chose to fix ^'^■^^(0) = ^l, 
j = 0, . . . , k — 2, which reduces in this case to fixing f2(0). We might, however, also want to 
optimize over this initial velocity. Unfortunately, examples can be exhibited where there does not 
exist any solution to the minimization of E if one also minimizes over f2(0). One possibility to 
restore well-posedness while retaining the minimization over $7(0) is to modify E by adding a 
penalty on the norm ||f2(0)||. 

The situation in infinite dimensions is similar, however proving existence results would require 
much deeper analytical study than in the finite dimensional case. 

In this section we presented higher-order methods that increase the smoothness in interpolating 
through a sequence of data points. In future work these methods will be compared to the shape 
spline model introduced in [TV 10]. Also of interest for CA is the metamorphosis approach that is 
discussed briefiy in Section 6. 



6 Optimization with penalty 

This section adapts the optimization approach of [GBHRIO] to higher-order Lagrangians. As in the 
case of the Clebsch-Pontryagin approach, one considers the (right or left) action ^ : G x Q ^ Q oi 
a Lie group G on a manifold Q. The basic idea is to replace the constraints in the Clebsch optimal 
control problem with a penalty function added to the cost function and to obtain in this way a 
classical (unconstrained) optimization problem. The penalty term is expressed with the help of a 
Riemannian metric 7 on the manifold Q. Given a cost function i : kg x Q ^ M., a > and the 
elements no, Ut € Q, C.^, E Q, j = 0, ...,k — 2, one minimizes 

£ (e, e, • • • , n) + ^Jn- ^Qin) f'^ dt, (6.1) 

over curves t ^ n(t) E Q and t ^(t) G 5 such that 

n(0)=no, n(T) = nT, ^^'\0) = ^^'\T) = ^t, J = 0,...,fc-2, 

where || ■ || is the norm on TQ induced by the metric 7 and, as in the Clebsch-Pontryagin case in 
Section 4, C,Q{n) denotes the infinitesimal generator of the G-action associated to ^ G 3, evaluated 
at n G Q. The corresponding stationarity conditions are found to be: 

k—l 

3=0 ^ ^ 



where the notation 



has been used and the covariant derivatives D/ Dt and V are associated to the Riemannian metric 
7 on Q. 

These equations should be compared with the stationarity conditions (4.2) associated to the 
Clebsch approach, which can be rewritten, with the help of a Riemannian metric, as 

g(-l)^9/^ = J(a), q = ^Q{q\ ^« = V^q) + |. (6.3) 
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Before proceeding further, we will pause to define some additional notation that will be conve- 
nient later. 

Definition 6.1. Consider a Lie group G acting on a Riemannian manifold {Q,'~f). We define the 
Q* -valued (1, 1) tensor field J-"^ : T*Q x TQ — q* associated to the Levi-Civita connection V hy 

{J^^iag, Ug),r]) := (a,, V„,r/Q(g)) , (6.4) 

for all Ug E TgQ, ag E T*Q, and rj E Q. 

The main properties of the tensor field J-"^ are discussed in [GBHRIO], where one also finds the 
proofs of the following two lemmas about the properties of J-"^. The first lemma below relates J-"^ 
to the connectors of the covariant derivatives on TQ and T*Q. The second lemma explains that 
J-"^ is antisymmetric under transposition in the inner product defined by the Riemannian metric 
7 when G acts by isometries. 

Lemma 6.2. For all ag E T*Q, Ug E TgQ, and ^ E Q, 

{j^^{ag,Ug),^) = {ag,K{^TQ{ug))) = - {K{^T*Q{aq)),Ug) , 
where K denotes the connectors of the covariant derivatives on TQ and T*Q, respectively. 

Proof. See the proof of Lemma 3.5 in [GBHRIO, §3]. □ 

Remark 6.3. A detailed treatment of connectors and their associated linear connections for co- 
variant derivatives can be found in [Mic08, §13.8]. We also refer to [GBHRIO] for useful properties 
of the connector K of relevance to the present paper. In infinite dimensions one needs to as- 
sume that the given weak Riemannian metric has a smooth geodesic spray 5* E X{TG), but such 
analytical issues will not be of concern to us here. 

Lemma 6.4. // G acts by isometries, then J-"^ is antisymmetric, that is 

J'^{ag,Uq) = -J'^(Mg,a5), 

for all Ug E TgQ, ag E T*Q. 

Proof. Since G acts by isometries, £(^Qg = 0; which implies (V.^q)^ = — V^q. □ 

The tensor field J-"^ arises naturally in computing the equations of motion associated to the 
stationarity conditions (6.2) for optimization with penalty. A computation, similar to the one 
given in [GBHRIO, §3] in the first order case, yields 

'"'^ „■ Si ^ fdi\ 1 



i=o ^ ^ ^ (6.5) 

D di 

where in (=l=) one chooses — (resp. +) when G acts on Q by a right (resp. left) action, consistently 
with (4.5). 

As a consequence of Lemma 6.4, if G acts on (Q, 7) by isometries, then the term ■^J-'^(u^, z/„) 
vanishes so that the optimization problem (6.1) produces the k'^^ order Euler-Poincare equations. 
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Example: The NHP equation via optimization. In this case, since S0{3) acts by isometrics 
on M.^, the minimization problem 



to 



min / {^\\nf + —\\q-nxqf]dt 



produces the NHP equations (1.1). 

Metamorphosis and Lagrange-Poincare reductions. Equations (6.5) may also be obtained 
by a generalization of the metamorphosis reduction developed in [GBHRIO], as follows. For sim- 
plicity, we only treat the case of a right action of G on Q. 

Consider a G-invariant Lagrangian L = L [g, g, g^''\ q, q) : T^'^^G x TQ ^ R relative to the 
action of /i G G given by 

{g,g,...,g^''\q,q) ^ {hg,hg, ...,hg^''\qh-\qh-^) 
and consider the quotient map 

{g, g, q, q) ^ (e, t ly) e kg x TQ, ^ = g-'g, n = qg,u = qg. (6.6) 

The equations of motion for the reduced Lagrangian ^M induced by L onkQX TQ can be obtained 
by a direct generalization of the method used in [GBHRIO] for k = 1. If L has the particular form 

L{g,g,...,g^''\q,q) = C {g, g, g'^''\q) + ^\\qgf, 

where C is the G-invariant Lagrangian associated to the function i in (6.1), then we recover 
equations (6.5) (with the upper sign chosen). 

Instead of the so-called metamorphosis quotient map (6.6) one may also use Lagrange-Poincare 
reduction with the quotient map 

{g,g,...,g^'\q,q)^ (^^,i...,C^'-'\n,f?j EkgxTQ, C = g-'g, n = qg. (6.7) 

The reduced equations of motion for metamorphosis with geometric splines that arise in the 
Lagrange-Poincare approach are 

' 51lp d6iLP 



6n dt 6h ' 

k~i (6.8) 
(a.TadJ)^(-iy9?|5^ = 0, 

j=0 ^ 

(with the upper sign chosen) where i^p is the reduced Lagrangian associated to the same unreduced 
Lagrangian L as before, but using the quotient map (6.7) instead of (6.6). 

Note that the Lagrange-Poincare approach generalizes easily to higher-order Lagrangians in q 
such as L := L [g,g, g^''\ q, q, : T'^^^G x Q) — t- M. The equations of motions are then 

simply 

i=o 

Si, 



j=0 
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The metamorphosis reduction approach also generahzes to higher higher-order Lagrangians in 
q. In this case, one uses the quotient map 

{g,g,...,g^'\q,q,...,q^'^) ^ (^^,t...,^^'-'\n,u,,...,u,^ekgxT^'^Q, (6.10) 

where ^ = g~^g and (ra, z/i, z/^) = ^g''^ (g, g, g'^'^^) , ^^'^^ being the natural induced action of 
G on T^^^Q. However for k > 2, the associated reduced equations are quite complex on general 
Riemannian manifolds so one may prefer to use the equivalent Lagrange-Poincare formulation 
(6.9). 

Remark 6.5. The idea of metamorphosis with splines may apply in imaging as in [HTY09] by 
using, e.g., L{gt, gt, guVt^Vt), L{gt, gt,VuVt, ,Vt), or L{gt, gt, gt,VuVt,Vt), instead of L{gt, gt^VuVt)- 



7 Clebsch and Lie-Poisson-Ostrogradsky formulations 

In this Section we present two Hamiltonian formulations associated to the higher order Euler- 
Poincare equations (3.29) with g-dependence. (The case of equations (3.30) may be obtained by 
making obvious modifications.) The first is a canonical Hamiltonian formulation that generalizes 
to higher order the canonical Clebsch formulation of Euler-Poincare dynamics. The second is 
a generalization of the Lie-Poisson formulation (with g-dependence) to higher order, that uses 
Ostrogradsky momenta. We now recall these formulations in the first order case. 

Clebsch canonical formulation. This is associated to the optimal control formulation de- 
scribed in §4. In the case k = 1 the canonical Hamiltonian formulation is already given by the 
Pontryagin approach. Indeed, if ^ || is a diffeomorphism we consider the function /i : g x Q — M 
defined by 

6i 

h{fi,q) := {^1,0 - T7=^ 

and the collective Hamiltonian H : T*Q — t- M given by H{aq) := h{3{aq)^ q). If aq{t) is a solution 
of Hamilton's canonical equations for H on T*Q, then {fi{t) , q{t)) , where fi{t) := J{ag{t)), is a 
solution of the Euler-Poincare equations 

This canonical formulation of the Euler-Poincare equations recovers some important examples 
such as the Clebsch variables for the ideal fluid [MW83], singular solutions of the Camassa-Holm 
equations [HM04], and double bracket equations, as explained in [GBR09]. 

Lie-Poisson formulation. This is obtained by reduction of the Hamiltonian Hg^ : T*G M. 
associated to Lg^ by Legendre transformation. If L is G- invariant as a function defined on TG x Q, 
then H : T*G x Q — )■ M is G-invariant and therefore induces the Hamiltonian h given above. By 
Poisson reduction of the manifold T*G x Q, where Q is endowed with the trivial Poisson structure, 
one obtains the Lie-Poisson equations 

(9,±ad|)/. = -j(^^), d,q-(^^^^{q) = 0, 
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together with the associated Poisson structure 

on g* X Q; see [GBTIO]. These equations are equivalent to their Lagrangian counterpart (7.1). 

7.1 Higher order Clebsch formulations 

Second order case 

Recall from Defintion 4. 1 that the Clebsch-Pontryagin variational formulation of the second order 
Euler-Poincare equations reads 



Sj'' {i{^,Lq) + {a,q-^Q{q)))dt = 0, 



over curves ^(t) G Q and a{t) G T*^^^Q and under conditions (A), (B), (C). If ^ i— t- tt := S£/6^ is a 

diffeomorphism, we define h{S,i'i^iq) '■= (/^, ^ K^y^^l) ^^"^ Pontryagin variational principle 
may be written equivalently as 

\^{tt,0 -h{^,7r,q) + {a,q-^Q{q))Jdt = 0, where it: 



5^ 



over curves ^(t) G and a{t) G T*^^^Q. Equivalently, this can be reformulated as 

6 j ((tt, i)-h (e, vr, q) + («, g - ^0(9))) dt = 0, 



over curves ^(t) G 0, 7r(t) G 0*, and a{t) G T*(j_-^Q, where 7r(i(:) is now an independent curve. The 

relation tt = Si/ 6^ is recovered by variations of One observes that this is simply the usual 

Hamilton Phase Space Variational Principle (i.e., not a Pontryagin Maximum Principle) on the 
phase space T*{Q x q) 

5 Hia,^,n)-(^{a,7r),{q,i)^dt = 0, 

for the Hamiltonian 

H:T*{Qxq)^R, H{a,, ^, tt) := h{^, tt, q) + (J(a,), ■ 
We thus have proved the following result. 

Theorem 7.1. Let £ : 20 x Q — ?■ M, £ = ^, q) be a cost function such that ,^ 1— t- tt := 5i/5^ is a 
diffeomorphism and define the function 

h{C,7C,q):={n,0-i{^,iq). (7.3) 

Then the stationarity conditions (4.4) for the 2"^^ order Clebsch-Pontryagin optimal control problem 
(4.1) with cost function i are given by the canonical Hamilton equations on T*{Q x 0) relative to 
the Hamiltonian H{^,fj,,a) = h{^,7i,q) + (J(a),^). 
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One can alternatively prove this result by computing explicitly the canonical Hamilton equa- 
tions for H on T*{Q x g). We obtain 

a = XH{a)=^T*Qia)-Vei^—, ^ = — = — , n = -— = --3{a). (7.4) 

Clearly, these equations coincide with the stationarity conditions (4.2). In particular, the last 
equation reads 

Example 1: Geodesic 2-spline equation on Lie groups 

Recall from §3.2 that the reduced Lagrangian for 2-splines on a Lie group G with right G-invariant 
Riemannian metric reads 

m,i) = l\\vr = l\\i + ^dl^f- (7.5) 

Here we denote 

7] = i + adje with adl V = (ad^(z/^))" for z/ G 0. (7.6) 
Then the quantity computed in equation (3.20) 



/i := ^ - dt-^ = [dtT] + ad^^ + adj^) , with r/ = ^ + adj^. 



satisfies the 2"''^-order Euler-Poincare equation, {dt + ad*^)fi = 0, which is also the geometric 2-spline 
equation of [CS95]. We now consider the canonical formulation of 2-splines. 

Hamiltonian formulation of the geodesic 2-spline equation on T*{Q x g). As we have 
seen, the Clebsch-Pontryagin approach of Section 4 allows the geodesic 2-spline equation to be 
recast as a set of canonical Hamilton equations for a Hamiltonian H : T*{Q x q) ^ M.. Note that 
in the case of 2-splines, the variable vr is 



5i 
5? 



vr = ^ = + ad*e^ = 7]^ 



which proves that ^ is a diffeomorphism. One thus obtains the Hamiltonian 

6i 

H{a, vr) = (tt, - K^, + (J(«), , ^ = 77 

= ^||7r|p-(7r,adJe) + (J(a),0, (7.7) 

where || ■ || denotes the norm induced by 7 on g*. The canonical Hamiltonian formulation (7.4) 
now yields the dynamical system 

« = ^T*Q{a), ^ = vr" - adj^, tt = -ad^j^^ - (ad^tt^^ - J («) (7.8) 

As we have proved above, the Euler-Poincare equation {dt + ad^/i = is then established by 
noticing that /i = J (a) is the cotangent-lift momentum map for the action of the Lie group G on 
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the manifold Q and that the ,^-equation impHes tt^ = ^ + ad^^ = rj. Note that the solution for 
the momentum map fi = J (a) may be obtained entirely from the canonical Hamilton equations, 
without explicitly solving the Euler-Poincare equation. For a bi-invariant metric, one has ad^,^ = 

in the ^-equation and the last two terms cancel each other in the 7r-equation. Consequently, these 
two canonical equations simplify to ^ = vr'* and tt = — J (a). From them, we find 

^ = -J(a)tt and = -adj^', (7.9) 

in agreement with equation (3.21) and reference [CS95]. 



Example 2: Geodesic 2-spline equations on 50(3) 

We consider the particular case of the Lie group G = S0{3) endowed with the bi-invariant metric 
induced by the standard Ad-invariant inner product 

7(l],r) = -^Tr(fir). 

We identify the dual 5o(3)* with 5o(3) using 7 so that fl^ = Q. Using the hat map ^: so(3) — )■ 
(see (3.23)), the Euler-Poincare equation in (7.9) reads 

h - n X n = 0, (7.10) 

which was first found in [NHP89]. The difference in sign from that paper arises here from the 
choice of reduction by right- invariance instead of left-invariance. 



Canonical Hamilton equations on T*]R^ x T*]R^ for the NHP equation. The Hamiltonian 
formulation of the NHP equation (7.10) for geometric splines on 5*0(3) with a bi-invariant metric 
may be obtained in canonical variables (f2,7r,q, p) G T*M.^ x T*M.^ from the Hamiltonian (see 
(7.7)), ^ 

i7(f2,7r,q,p) = ^||7rf + ri-qxp. (7.11) 

This corresponds to the choice Q = M'^ on which 50(3) acts by matrix multiplication. 
This Hamiltonian produces canonical equations of the form, 

q= =f2xq, p = -— = rixp, 
op oq 

^ SH . SH 

n= — =7Z, 77 = -— = -qxp. 

OTT Oil 

The (q, p)-equations here imply that = J(q, p) = q x p obeys the Euler-Poincare equation for 
right invariance, 

{dt + ad^) n = = fi — Qxfi, 

which results in the NHP equation (7.10) when we substitute ^ = —Cl. 

The canonical Hamiltonian formulation of the NHP equation provides some insight into the 
interpretation of its constants of motion. For example, the Hamiltonian (7.11) Poisson commutes 
with |qp, IpP, (q ■ p), and |q x pp, although only the last of these Poisson commutes with all the 
others. The Hamiltonian (7.11) also Poisson commutes with the vector K = f2x7r — qxpG 
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M^. Although the components of K satisfy Poisson bracket relations {Ki,K2} = and cyclic 
permutations with each other, their sum of squares = K\ + fr| + again Poisson commutes 
with all the others. The presence of the two constants of motion |q x pp and K"^ in Poisson 
involution allows symplectic reduction from six degrees of freedom to four, but the reduced system 
is still far from being integrable. The Hamiltonian conservation laws may be expressed in terms 
of (O, O, n) e T(2)m3 as 

|q X p|2 = K = rixr2 + fi, k'^ = \nxn\'^ + 2{nxn) ■n + ini'^. 

All of these conservation laws were known in the literature, but had previously not been given a 
Hamiltonian interpretation. The Hamiltonian interpretation of the NHP equation (7.10) in this 
setting is that the rotations act on the cross product m = q x p diagonally in q and p, so that 
m = f2xmforq = f2xq and p = 17 x p. This is also the essence of the symmetric representation 
of rigid body motion discussed, e.g., in [BC96b]. 

Higher order case 

The canonical Clebsch formulation presented above can be adapted to higher order cost functions 

^ jg(fc-i) is a diffeomorphism, we define the function 

■■= {r^2,0 + (^3,0 + - + (vrfc,^^'^-^)) - vr, = (7.12) 

and we consider the Hamiltonian H : T*{Q x {k — 1)q) — )■ M given by 

H{a„^,...,^^''-'\n2,...,7rk) := /i (e, e, vr^, vr,, g) + (J(a,), • 

A straightforward computation shows that the canonical Hamilton equations on T*{Q x {k — 1)q) 
for H produce the stationarity condition (7.1) of the fc*'*-order Clebsch- Pontryagin optimal control 
with cost function i and therefore imply the A;*^-order Euler-Poincare equations (3.29). We thus 
obtain the generalization of Theorem 7.1 for k^'^-oidei cost functions. 

7.2 Ostrogradsky-Lie-Poisson reduction 

The procedure of Lie-Poisson reduction of the Hamilton-Ostrogradsky theory parallels that for 
higher-order Euler-Poincare reduction and produces a different Hamiltonian formulation of the 
higher-order dynamics that applies to A; > 2. At first, we will discuss the Hamilton-Ostrogradsky 
approach for the higher-order Hamiltonian formulation based purely on Lie group reduction, i.e., 
without introducing the action of the Lie group G on the manifold Q. Then we will remark on 
how g-dependence may be easily incorporated. 

Second order. Consider a G-invariant second order Lagrangian L : T^'^^'G — t- M, L = L{g,g,g). 
The Ostrogradsky momenta are defined by the fiber derivatives. 



dL ^ dL dL 

Pl = Ot — , P2 = ^ 

og og og 
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and produce the Legendre transform {g,g,g, 9) E T^^'^G t— ?■ {g, g,Pi,P2) ^ T*(TG). We refer to 
[dLR85] for the intrinsic definition of the Legendre transform for higher order Lagrangians. See 
also [BC96a] for an apphcation on 5*0(3). 

When the Legendre transform is a diffeomorphism, the corresponding Hamiltonian H : T*{TG) — )■ 
M is defined by 

H{g,g,Pi,P2) ■■= {pug) + {P2,9) - L{g,g,g) 

and the canonical Hamilton equations for H are equivalent to the 2"''^-order Euler-Lagrange equa- 
tions for L. 

Applying reduction by symmetry to H induces a Hamiltonian /i(7ri, 7^2) on T*{TG) /G ~ g* x 
T*g, which is related to the symmetry-reduced Lagrangian £(^, S,) by the corresponding Legendre 
transformation, 



/i(7ri, e, vrs) = (vri, + (tts, ^ - £(e, 0, ^ = ^2- (7.13) 

By Reduction of the Hamilton-Ostrogradsky equations for H on T* {TG) we obtain the Ostrogradsky- 
Lie-Poisson equations for h 



dtTTi ± ad*^fe TTi = 0, 

n (7.14) 

6712 



together with the non-canonical Poisson bracket given by 
{/,5'}(7ri,^,vr2) = ± ( VTi, 



5/ 5/ 
{/,^7}±(vri) + {/,^/}can(e,vr2) (7.15) 



6ni ' ^TTi 



for functions f,g depending on the variables (7ri,^,7r2). Note that this reduction process holds 
without assuming a preexisting Lagrangian formulation. Equations (7.14) together with their 
Hamiltonian structure can be obtained by Poisson reduction for cotangent bundles: T*Q — )■ T*Q/G 
(the so called Hamilton-Poincare reduction [CMPR03]) applied here to the special case Q = TG. 

We now check directly that equations (7.14) are equivalent to the 2"''^-order Euler-Poincare 
equations if the Hamiltonian (7.13) is associated to i by an invertible Legendre transform. The 
derivatives of the symmetry-reduced Hamiltonian h with respect to the momenta vri and 7^2 imply 
formulas for the velocity and acceleration, 

^ Sh ■ X 

so that the acceleration ^ may be expressed as a function of the velocity ^ and the momenta 
(7ri,7r2). The pair (^, vr2) G T*g ^ g x g* obeys canonical Hamilton equations, so the derivatives 
of h in (7.13) with respect to velocity and acceleration imply the momentum relations. 



(7.17) 
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Solving these momentum relations for tti and tt2 in terms of derivatives of the reduced Lagrangian 
yields 

TTi = — - dt-j and TTa = -j. (7.18 

The Lie-Poisson equation for tti, 

StTTi ± ad*«^ TTi = 0, (7.19) 
then implies the 2"'^-order Euler-Poincare equation, 

Example: Ostrogradsky-Lie-Poisson approach for geometric 2-splines 

The Ostrogradsky reduced Hamiltonian (7.13) for geometric 2-splines is 

/i(7ri,e,7r2) = ^ ||7r2f - (7r2,ad^e) + (^i,^)- (7.20) 

From this reduced Hamiltonian, the Poisson bracket (7.15) recovers the geometric 2-spline equa- 
tions (3.20). For a bi-invariant metric ad^.^ = and these equations reduce to (3.21). In addition 
for 50(3) these equations produce the NHP equation (7.10). 

Third order. Before going to the general case, it is instructive to quickly presenting the case 
of a third order G-invariant Lagrangian L : T^'^^G — )■ M, L = L(g,g,g, 9) inducing the symmetry- 
reduced Lagrangian i : T^^^G/G ~ 30 ^ M, £ = 1 0- The Ostrogradsky momenta 

dL ^ dL ^2 dL dL ^ dL dL 

og og dg og dg dg 

produce the Legendre transform (g, g, g^^^) E T^^^G t-)- {g, g, 'g,Pi,P2,P3) £ T*(T^'^^G). The 
associated G-invariant Hamiltonian is obtained from the Legendre transformation, 

H : r*(T(^)G) ^ M, Hig,g,g,p,,p2,P3) ■= {pi,g) + {p2,g) + {ps, 'd) - Lig,g,g,g), 

so that G-invariance of the Hamiltonian H yields the symmetry-reduced Hamiltonian h{TTi, ^, ^, 7r2, tts), 
h : T*{T^'^'>G)/G ~ 0* X T*(20) M. The reduced Hamiltonian h is related to the reduced La- 
grangian ^ by the extended Legendre transformation 

/i(7ri, e, L 7r2, TTs) = (tti, + (vra, i) + (vra, 'C) - ^{i, 1 0, ^ = ^3- (7.21) 
The 3'''^-order Ostrogradsky-Lie-Poisson system reads 



( dtiii ± ad*«h TTi = 0, 


















(7.22) 
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with associated Poisson bracket 

{/, gU^ru t vr2, TTa) = {/, h}i{n,) + {/, n^) + {/, hj^^^it tt,). (7.23) 

If the Hamiltonian (7.21) is associated to a Lagrangian i by Legendre transformation, we have 
■ 6h . 6h 6£ ■■ 6h . 6h 6£ . ^ 

4 = —, VTs = -— = -7ri + — , 4 = ^, TTg = -— r = -7r2 + -7, TTi + ad^ TTi = 0. (7.24) 
dTT2 Ot, Ot, dTTs 6^ 

Consequently, we have 

and the last equation in (7.24) for the momentum map vri imphes by the 3rd-order Euler-Poincare 
equation 



Higher-order and g-dependence. The Ostrogradsky-Lie-Poisson approach generalizes to k^'^- 
order as follows. For a G-invariant Lagrangian L : T^'^^G — )■ M, L = L(g, g, g^'^^), the Ostrograd- 
sky momenta define the Legendre transform as a map T^'^''~^^G — > T*{T^'^~^^G) (see [dLR85]) and 
the associated Hamiltonian H : T*{T^'^~^^G) — )■ M, if = H(g,g, g^'^^^\pi, ...,Pk) is given by 

k 

H{g,g, ...,g'^'''\pu . . . ,p,) := ^ {g(^\p,) - L{g,g, ...,g^'^). (7.25) 

i=i 

Extending the symmetry-reduced Ostrogradsky procedure outlined above to /c^^'-order then leads 
to the Ostrogradsky-Lie-Poisson equations 




whose Hamiltonian structure is 



k 



{/, ^?}(7ri, . . . , vrs, Hk) = {/, ^7}±(7ri) + g].an{i^^~'\ vr,). (7.27) 

i=2 

This Poisson bracket produces the geometric fc-spline equations from the corresponding reduced 
Hamiltonian; we do not carry out the details here. 

The Ostrogradsky procedure for reduction by symmetry outlined above generalizes to allow 
g-dependence. Indeed, for a G-invariant Lagrangian L : {T^^'^G) X g ^ M, the previous steps may 
all be repeated with only slight changes, resulting in the reduced Poisson bracket (7.27), modified 
by adding the terms 
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Thus, allowing g-dependence leaves the canonical equations invariant, but alters the vri-equation 
so that it becomes 

dtTTi ± adV TTi = - J ( ^ ) . (7.29) 



Prom Clebsch to Ostrogradsky-Lie-Poisson. Recall that, in the case k = 1, one passes from 
the canonical Clebsch formulation to the Lie-Poisson formulation (with Lie-Poisson bracket (7.2)) 
by using the momentum map, via the transformation 

ag G T*Q (/i, q) eg* xQ, jj = 3{aq). 

Consider now the case k = 2 with Lagrangian i = ^, q). On the Clebsch side, the Hamiltonian 
is given by 



H{ag, e, tt) := (vr, - i q) + (J(«,), , 




whereas the Ostrogradsky-Lie-Poisson Hamiltonian with g-dependence is defined by 

/i(7ri, ^, VTa, q) := (tti, ^) + {^773, t) - ^(^, ^, ?), = ^ • 

These definitions suggest that one can pass from the canonical Clebsch formulation on T*{Q x q) 
to the Lie-Poisson-Ostrogradsky formulation on g* x T*q x Q by the Poisson map 

(ttg, ^, TTa) (tti, ^, TTa, g), TTi := J (a,). 

This is indeed the case, as one can check easily. Generalization to A; > 2 is now straightforward 
and we have the Poisson map 

T*{Q x{k- l)s) ^ s* X T*{k - 1)0 X Q, 

that relates the canonical Hamilton equations on T*{Qx (A; — 1)^) and the Ostrogradsky-Lie-Poisson 
equations (7.26). 



8 Outlook and open problems 

8.1 Brief summary and other potential directions 

This paper has begun the application of symmetry-reduction tools to higher-order variational 
problems on Lie groups, culminating in an application of 2"'^-order geometric splines to template 
matching on the sphere which was shown to be governed by a higher-order Euler-Poincare equation 
on the dual Lie algebra of the Lie group SO (3). The generality of this result was emphasized in 
Remark 5.2, on seeing that the higher-order Euler-Poincare equation (5.15) had emerged once 
again as the optimality condition for template matching. 

Various open problems not treated here seem to crowd together to present themselves. A few 
of these are: 
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• We have applied variational constraints to /c-sp lines in Section 3.4. However, accommodating 
nonholonomic constraints would require additional developments of the theory. 

• In Section 5 we presented higher-order methods that increased the smoothness in interpo- 
lating through a sequence of data points. In future work these methods will be compared to 
the shape spline model introduced in [TVIO]. Some initial forays into the analysis of these 
problems were also presented in Section 5.2, but much remains to be done for these problems 
that have been treated here only formally. 

• Extension of the basic theory presented here to allow for actions of Lie groups on Rieman- 
nian manifolds is also expected to have several interesting applications, particularly in image 
registration. For example, one could address higher-order Lie group invariant variational 
principles that include both curves on Lie groups and the actions of Lie groups on smooth 
manifolds, particularly on Riemannian manifolds. These Lie group actions on manifolds ap- 
ply directly to the optimal control problems associated with large-deformation image regis- 
tration in the Large Deformation by Diffeomorphisms Metric Mapping (LDDMM) framework 
via Pontryagin's maximum principle. Actions of Lie groups on Riemannian manifolds will 
be investigated in a subsequent treatment. 

Doubtlessly, other opportunities for applying and extending this symmetry reduction approach 
for /c-splines will present themselves in further applications. 

8.2 An open problem: the slalom, or brachistochrone for splines 

Let us formulate yet another example in slightly more detail. This is the brachistochrone version 
of the optimization problem treated here, for possible application, say, in a slalom race. Unlike 
the optimization problem, which seeks the path of least cost, a race would seek the path of least 
time. For example, the familiar slalom race involves dodging around a series of obstacles laid 
out on the course. The objective of the slalom racer in down-hill skiing, for example, is to pass 
through a series of gates as quickly as possible. The strategy in slalom racing is to stay close to 
the shortest-time path (or geodesic) between the gates, while also moderating the force exerted in 
turning to keep it below some threshold, lest the snow give way and the skier slides off the course. 
Thus, the ideal slalom path sought by an expert racer would hug the geodesic between the gates 
and make the series of turns passing through the gates with no skidding at all. 

The strategy for achieving the optimal slalom has many potential applications in modern 
technology. For example: 

• A charged-particle beam in an accelerator may be guided in its path by a series of quadrupole 
magnets that steer the beam to its target. The steering must be done as gently as possible, 
so as to minimize the transverse acceleration (seen as curvature in the path of the beam) 
that causes Bremsstrahlung and the consequent loss of energy in the beam. 

• An underwater vehicle may be steered smoothly through narrow passageways in a sunken 
ship along a path that will take it quickly and efficiently to its objective, thereby avoiding 
collisions while minimizing fuel expenditure, time, etc. 

• A car may be programmed to glide smoothly through a tight parallel-parking maneuver that 
ends in an elegant stop in the narrow space between two other cars along the curb. 
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• A vehicle may follow a program to roll as rapidly as possible along the terrain through a 
series of gates with its cameras mounted so that they continuously point toward an object 
above it that must keep in sight. 

The slalom strategy that apphes in all these examples seeks a path that minimizes the time 
for the distance travelled over a prescribed course, while also moderating the acceleration or force 
exerted along the path as it passes around a series of obstacles or through a series of gates laid 
out on the course. Designing such maneuvers requires optimization for least time, while also using 
cost functions that depend on both the velocity and acceleration of the motion. Moderation of 
higher-order accelerations such as jerk (rate of change of acceleration) may also be needed. As in 
the present paper, in solving optimal slalom problems that minimize the time taken to finish the 
course, one might expect to take advantage of continuous symmetries by investigating Lie group 
invariant variational problems for cost functions that are defined on A;*^-order tangent spaces of Lie 
groups acting on smooth Riemannian manifolds. Investigations of invariant variational principles 
for slalom problems using the present group reduction and induced metric methods would be 
a promising direction. However, this direction seems even more challenging than the geometric 
splines for optimizing costs in trajectory planning on a Lie group treated here and it will be deferred 
to a later paper. 
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